﻿//TODO: test, Gabor

using System;
using System.Drawing;
using System.Drawing.Drawing2D;
using System.Drawing.Imaging;
using FourierTools.DSP;
namespace LowLevelGraphics.Filter
{
    /// <summary>
    /// Gabor Transformation
    /// </summary>
    public class Gabor : BaseImageFilter
    {
        private int[,] m_aImage = null;

        protected int m_nOrientation = 0;
        protected int m_nFrequency = 0;
        protected int m_nWindowSize = 256;

        protected int m_nGaborWidth = 101;
        protected int m_nGaborHeight = 101;
        private Complex[] KernelFFT2 = null;

        /// <summary>
        /// Initializes a new instance of the <see cref="Gabor"/> class.
        /// </summary>
        public Gabor()
        {
            InitializeVariables();
        }

        /// <summary>
        /// Initialize Gabor Transfomation Window
        /// </summary>
        private void InitializeVariables()
        {
            KernelFFT2 = new Complex[m_nWindowSize * m_nWindowSize];
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="Gabor"/> class.
        /// </summary>
        /// <param name="_gabor">The _gabor.</param>
        internal Gabor(Gabor _gabor)
        {
            this.m_aImage = _gabor.m_aImage;
        }

        /// <summary>
        /// Gabor Width
        /// </summary>
        public int GaborWidth
        {
            get { return m_nGaborWidth; }
            set { m_nGaborWidth = value; }
        }

        /// <summary>
        /// Gabor Height
        /// </summary>
        public int GaborHeight
        {
            get { return m_nGaborHeight; }
            set { m_nGaborWidth = value; }
        }

        /// <summary>
        /// Sets or gets the Frequency
        /// </summary>
        public int Frequency
        {
            get { return m_nFrequency; }
            set { m_nFrequency = value; }
        }

        /// <summary>
        /// Sets or gets the orientation
        /// </summary>
        public int Orientation
        {
            get { return m_nOrientation; }
            set { m_nOrientation = value; }
        }

        /// <summary>
        /// Gets or sets Window Size
        /// </summary>
        public int WindowSize
        {
            get { return m_nWindowSize; }
            set { m_nWindowSize = value; }
        }

        /// <summary>
        /// Execute this Filter
        /// A Debug Statement is written if a concrete Filter is called
        /// </summary>
        /// <param name="_bitmap"></param>
        /// <returns></returns>
        public override UnsafeBitmap Execute(UnsafeBitmap _bitmap)
        {
            UnsafeBitmap bitmap = _bitmap;
            PrepareImage(bitmap);
            PrepareKernel();
            GaborResult result = new GaborResult();
            GaborTransform(_bitmap, m_aImage, m_nOrientation, m_nFrequency, result);
            return _bitmap;
        }

        /// <summary>
        /// Prepares the image.
        /// </summary>
        /// <param name="_bitmap">The _bitmap.</param>
        protected void PrepareImage(UnsafeBitmap _bitmap)
        {
            int nHeight = _bitmap.Height;
            int nWidth = _bitmap.Width;

            m_aImage = new int[nHeight, nWidth];

            for (int y = 0; y < nHeight; y++)
            {
                for (int x = 0; x < nWidth; x++)
                {
                    m_aImage[y, x] = ColorToGray(_bitmap.GetPixel(x, y));
                }
            }
        }

        /// <summary>
        /// Colors to gray.
        /// </summary>
        /// <param name="color">The color.</param>
        /// <returns></returns>
        protected int ColorToGray(Color color)
        {
            //Gray = 0.299 * R + 0.587 * G + 0.114 * B
            //Gray = (R * 77 + G * 151 + B * 28) >> 8; 
            int gray = ((((color.R * 77) + (color.G * 151)) + (color.B * 28)) >> 8);

            if (gray < 0) gray = 0;
            if (gray > 255) gray = 255;

            return gray;
        }

        /// <summary>
        /// Prepares the kernel.
        /// </summary>
        public void PrepareKernel()
        {
            CalculateKernel(m_nOrientation, m_nFrequency);
            Fourier.FFT2(KernelFFT2, m_nWindowSize, m_nWindowSize, FourierDirection.Forward);
        }

        /// <summary>
        /// Calculates the kernel.
        /// </summary>
        /// <param name="_nOrientation">The _n orientation.</param>
        /// <param name="_nFrequency">The _n frequency.</param>
        private void CalculateKernel(int _nOrientation, int _nFrequency)
        {
            double real, img;

            for (int x = -(m_nGaborWidth - 1) / 2; x < (m_nGaborWidth - 1) / 2 + 1; x++)
            {
                for (int y = -(m_nGaborHeight - 1) / 2; y < (m_nGaborHeight - 1) / 2 + 1; y++)
                {
                    real = KernelRealPart(x, y, _nOrientation, _nFrequency);
                    img = KernelImgPart(x, y, _nOrientation, _nFrequency);

                    KernelFFT2[(x + (m_nGaborWidth - 1) / 2) + m_nWindowSize * (y + (m_nGaborHeight - 1) / 2)].Real = real;
                    KernelFFT2[(x + (m_nGaborWidth - 1) / 2) + m_nWindowSize * (y + (m_nGaborHeight - 1) / 2)].Imaginary = img;
                }
            }
        }

        /// <summary>
        /// Kernels the real part.
        /// </summary>
        /// <param name="x">The x.</param>
        /// <param name="y">The y.</param>
        /// <param name="_nOrientation">The _n orientation.</param>
        /// <param name="_nFrequency">The _n frequency.</param>
        /// <returns></returns>
        protected double KernelRealPart(int x, int y, int _nOrientation, int _nFrequency)
        {
            double U, V;
            double dSigma, Kv, Qu;
            double tmp1, tmp2;

            U = _nOrientation;
            V = _nFrequency;
            dSigma = 2 * Math.PI * Math.PI;
            Kv = Math.PI * Math.Exp((-(V + 2) / 2) * Math.Log(2, Math.E));
            Qu = U * Math.PI / 8;

            tmp1 = Math.Exp(-(Kv * Kv * (x * x + y * y) / (2 * dSigma)));
            tmp2 = Math.Cos(Kv * Math.Cos(Qu) * x + Kv * Math.Sin(Qu) * y) - Math.Exp(-(dSigma / 2));

            return tmp1 * tmp2 * Kv * Kv / dSigma;
        }

        /// <summary>
        /// Kernels the img part.
        /// </summary>
        /// <param name="x">The x.</param>
        /// <param name="y">The y.</param>
        /// <param name="_nOrientation">The _n orientation.</param>
        /// <param name="_nFrequency">The _n frequency.</param>
        /// <returns></returns>
        protected double KernelImgPart(int x, int y, int _nOrientation, int _nFrequency)
        {
            double U, V;
            double Sigma, Kv, Qu;
            double tmp1, tmp2;

            U = _nOrientation;
            V = _nFrequency;
            Sigma = 2 * Math.PI * Math.PI;
            Kv = Math.PI * Math.Exp((-(V + 2) / 2) * Math.Log(2, Math.E));
            Qu = U * Math.PI / 8;

            tmp1 = Math.Exp(-(Kv * Kv * (x * x + y * y) / (2 * Sigma)));
            tmp2 = Math.Sin(Kv * Math.Cos(Qu) * x + Kv * Math.Sin(Qu) * y) - Math.Exp(-(Sigma / 2));

            return tmp1 * tmp2 * Kv * Kv / Sigma;
        }

        //TODO: interpret results back to IMAGE
        /// <summary>
        /// Gabors the transform.
        /// </summary>
        /// <param name="_bitmap">The _bitmap.</param>
        /// <param name="_aImage">The _a image.</param>
        /// <param name="_nOrientation">The _n orientation.</param>
        /// <param name="_nFrequency">The _n frequency.</param>
        /// <param name="_result">The _result.</param>
        public void GaborTransform(UnsafeBitmap _bitmap, int[,] _aImage, int _nOrientation, int _nFrequency, GaborResult _result)
        {
            Complex[] imageFFT = new Complex[m_nWindowSize * m_nWindowSize];
            Complex[] ConvolutionResult = new Complex[m_nWindowSize * m_nWindowSize];

            int nWindowSizeHalf = m_nWindowSize / 2;
            for (int y = (m_nGaborWidth - 1) / 2; y < (m_nGaborWidth - 1) / 2 + nWindowSizeHalf; y++) //            for (int y = (GaborWidth - 1) / 2; y < (GaborWidth - 1) / 2 + 128; y++)
            {
                for (int x = (m_nGaborWidth - 1) / 2; x < (m_nGaborWidth - 1) / 2 + nWindowSizeHalf; x++) //for (int x = (GaborWidth - 1) / 2; x < (GaborWidth - 1) / 2 + 128; x++)
                {
                    imageFFT[y * m_nWindowSize + x].Real = _aImage[y - (m_nGaborWidth - 1) / 2, x - (m_nGaborWidth - 1) / 2];
                }
            }

            Fourier.FFT2(imageFFT, m_nWindowSize, m_nWindowSize, FourierDirection.Forward);

            for (int i = 0; i < m_nWindowSize * m_nWindowSize; i++)
                ConvolutionResult[i] = imageFFT[i] * KernelFFT2[i] / (m_nWindowSize * m_nWindowSize);

            Fourier.FFT2(ConvolutionResult, m_nWindowSize, m_nWindowSize, FourierDirection.Backward);

            double dSum = 0, dAvg = 0, dDeta = 0;
            double dTempModulus = 0;
            double[,] aTempMag = new double[m_nWindowSize, m_nWindowSize];

            for (int y = (m_nGaborWidth - 1); y < (m_nGaborWidth - 1) + nWindowSizeHalf; y++) //for (int y = (GaborWidth - 1); y < (GaborWidth - 1) + 128; y++)
            {
                for (int x = (m_nGaborWidth - 1); x < (m_nGaborWidth - 1) + nWindowSizeHalf; x++) // for (int x = (GaborWidth - 1); x < (GaborWidth - 1) + 128; x++)
                {
                    dTempModulus = ConvolutionResult[y * m_nWindowSize + x].GetModulus();

                    dSum += dTempModulus;
                    aTempMag[y, x] = dTempModulus;
                }
            }

            dAvg = dSum / (nWindowSizeHalf * nWindowSizeHalf); //  dAvg = dSum / (128 * 128);

            for (int y = (m_nGaborWidth - 1); y < (m_nGaborWidth - 1) + nWindowSizeHalf; y++) // for (int y = (GaborWidth - 1); y < (GaborWidth - 1) + 128; y++)
                for (int x = (m_nGaborWidth - 1); x < (m_nGaborWidth - 1) + nWindowSizeHalf; x++) //for (int x = (GaborWidth - 1); x < (GaborWidth - 1) + 128; x++)
                    dDeta += (aTempMag[y, x] - dAvg) * (aTempMag[y, x] - dAvg);

            dDeta = dDeta / (nWindowSizeHalf * nWindowSizeHalf); // dDeta = dDeta / (128 * 128);

            for (int y = (m_nGaborWidth - 1) + 0; y < (m_nGaborWidth - 1) + 8; y++)
            {
                for (int x = (m_nGaborWidth - 1) + 0; x < (m_nGaborWidth - 1) + 8; x++)
                {
                    Console.Write("{0,8:F2}", aTempMag[y, x]);
                }
                Console.WriteLine();
            }

            for (int y = 0; y < m_nWindowSize; y++)
            {
                for (int x = 0; x < m_nWindowSize; x++)
                {
                    int tmpMag1 = (int)aTempMag[y, x];
                    if (tmpMag1 > 255) tmpMag1 = 255;
                    _bitmap.SetPixel(x, y, Color.FromArgb(tmpMag1, tmpMag1, tmpMag1));
                }
            }
            _result.Avg = dAvg;
            _result.Deta = dDeta;
        }

        /// <summary>
        /// This is used to clone a filter, each filter must implement it's own clone functionality
        /// </summary>
        /// <returns></returns>
        public override object Clone()
        {
            return new Gabor(this);
        }
    }
}
